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This article proposes a topological method that extracts hierarchi¬ 
cal structures of various amorphous solids. The method is based 
on the persistence diagram (PD), a mathematical tool for captur¬ 
ing shapes of multi-scale data. The input to the PDs is given by 
an atomic configuration and the output is expressed as two dimen¬ 
sional histograms. Then, specific distributions such as curves and 
islands in the PDs identify meaningful shape characteristics of the 
atomic configuration. Although the method can be applied to a wide 
variety of disordered systems, it is applied here to silica glass, the 
Lennard-Jones system, and Cu-Zr metallic glass as standard exam¬ 
ples of continuous random network and random packing structures. 
In silica glass, the method classified the atomic rings as short-range 
and medium-range orders, and unveiled hierarchical ring structures 
among them. These detailed geometric characterizations clarified 
a real space origin of the first sharp diffraction peak, and also in¬ 
dicated that PDs contain information on elastic response. Even in 
the Lennard-Jones system and Cu-Zr metallic glass, the hierarchical 
structures in the atomic configurations were derived in a similar way 
using PDs, although the glass structures and properties substantially 
differ from silica glass. These results suggest that the PDs provide a 
unified method that extracts greater depth of geometric information 
in amorphous solids than conventional methods. 

Amorphous Solid | Hierarchical Structure | Persistence Diagram | Topological 
Data Analysis 

Abbreviations; PD, persistence diagram; SRO, short-range order; MRO, medium- 
range order; FSDP, first sharp diffraction peak 

T he atomic configurations of amorphous solids are diffi¬ 
cult to characterize. Because they have no periodicity as 
found in crystalline solids, only local structures have been an¬ 
alyzed in detail. Although short-range order (SRO) defined by 
the nearest neighbor is thoroughly studied, it is not sufficient 
to fully understand the atomic structures of amorphous solids. 
Therefore, medium-range order (MRO) has been discussed to 
properly characterize amorphous solids [1, 2, 3]. Many experi¬ 
mental and simulation studies [4, 5, 6, 7] have suggested signa¬ 
tures of MRO such as a first sharp diffraction peak (FSDP) in 
the structure factor of the continuous random network struc¬ 
ture, and a split second peak in the radial distribution func¬ 
tion of the random packing structure. However, in contrast to 
SRO, the geometric interpretation of MRO and the hierarchi¬ 
cal structures among different ranges are not yet clear. 

Among the available methods, the distributions of bond an¬ 
gle and dihedral angle are often used to identify the geometry 
beyond the scale of SRO. They cannot, however, provide a 
complete description of MRO because they only deal with the 
atomic configuration up to the third nearest neighbors. Alter¬ 
natively, ring statistics are also applied as a conventional com¬ 
binatorial topological method [2, 8, 9]. However, this method 
is applicable only for the continuous random network or crys¬ 
talline structures, and furthermore, cannot describe length 
scale. Therefore, methodologies that precisely characterize hi¬ 
erarchical structures beyond SRO and are applicable to a wide 
variety of amorphous solids are highly desired. 

In recent years, topological data analysis [10, 11] has 
rapidly grown and has provided with several tools for study¬ 


ing multi-scale data arising in physical and biological fields 
[11, 12, 13, 14, 15, 16]. A particularly important tool in the 
topological data analysis is persistence diagram (PD), a vi¬ 
sualization of persistent homology as a two-dimensional his¬ 
togram (e.g., see Fig. 2). The input to the PD is given by an 
atomic configuration with scale parameters, and the output 
consists of various multi-scale information about topological 
features such as rings and cavities embedded in the atomic 
configuration. Here, the atomic configurations are generated 
by molecular dynamics simulations in this article. Impor¬ 
tantly, in contrast to other topological tools, PDs not only 
count topological features, but also provide the scales of these 
features. Hence, PDs can be used to classify topological fea¬ 
tures by their scales and clarify geometric relationships among 
them; this is presumably the most desired function for deeper 
analysis of amorphous structures. 

This article proposes a method using PDs for various amor¬ 
phous solids in a unified framework. The method is applied 
to atomic configurations and enables one to study hierarchi¬ 
cal geometry embedded in amorphous structures, which can¬ 
not be treated by conventional methods. We first applied the 
method to silica glass as an example of the continuous ran¬ 
dom network structure, and obtained the following results: (1) 
We found three characteristic curves in the PD of silica glass. 
These curves classify the SRO rings in the Si 04 tetrahedra 
and the MRO rings constructed by those tetrahedra. Further¬ 
more, a hierarchical relationship among the SRO and MRO 
rings was elucidated. (2) The PD reproduced the wavelength 
of the FSDP and clarified a real space origin of the FSDP. 
(3) Each curve in the PD represents a geometric constraint on 
the ring shapes and, as an example, an MRO constraint on 
rings consisting of three oxygen atoms was explicitly derived 


Significance 

Persistent whomology is an emerging mathematical concept for 
characterizing shapes of data. In particular, it provides a tool 
called the persistence diagram that extracts multi-scale topo¬ 
logical features such as rings and cavities embedded in atomic 
configurations. This article presents a unified method using per¬ 
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rations in amorphous solids. The method highlights hierarchical 
structures that conventional techniques could not have treated 
appropriately. 
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as a surface in a parameter space of the triangles. Moreover, 
we verified that these curves are preserved under strain, in¬ 
dicating that the PD properly encodes the material property 
of elastic response. Next, as examples of the random pack¬ 
ing structure, the Lennard-Jones system and Cu-Zr metallic 
glass were studied by the PDs, and we clarified the follow¬ 
ing: (4) These amorphous solids were also characterized well 
by the distributions of curves and islands in the PDs. (5) In 
the Lennard-Jones system, the global connectivity of dense 
packing regions was revealed by dualizing octahedral arrange¬ 
ments. (6) In Cu-Zr alloys, we found that the pair-distribution 
function defined by the octahedral region in the PD shows the 
split second peak. Furthermore, a relationship between the 
hierarchical ring structure and high glass forming ability was 
discovered in Cu-Zr alloys. 


PDs of atomic configuration 

The input to PDs is a pair A = (Q, B,) of an atomic configura¬ 
tion Q = (xi,... ,X]v) and a parameter set R = (^i, • • • 

Here, Xi and are the position in and the input ra¬ 
dius for the ith atom, respectively. In order to character¬ 
ize the multi-scale properties in Q, we introduce a param¬ 
eter a, which controls resolution, and generate a family of 
atomic balls Bi(a) = {xGR^ | ||cc — a3i|| < ri{a)} having the 
radius ri{a) = arf. We vary the radii of the atomic 
balls by a and detect rings and cavities at each a, where 
a > amin := - min{ri,..., r%}. 

Let Ck be a ring or cavity constructed in the atomic ball 
model B{a) = U^i ^ parameter a. To be more 

precise, a ring (resp. cavity) here means a generator of the 
homology Hi{B{a)) (resp. H 2 {B{a))) with a field coefficient 
[17]. Then, we observe that there is a value a — bk (resp. 
a = d/c) at which Ck first appears (resp. disappears) in the 
atomic ball model. The values bk and dk are called the birth 
and death scales of Cfc, respectively. The collection of all the 
(bk,dk) G R^ of rings (resp. cavities) is the PD denoted by 
Di{A) (resp. D 2 {A)) for A (see Fig. I). It follows from the 
structure theorem of persistent homology [II] that the PDs 
are uniquely defined from the input. From this construction, 
(bk,dk) encodes certain scales of each Ck- For example, in 
Di{A), bk indicates the maximum distance between two ad¬ 
jacent atoms in the ring c^, whereas dk indicates the size of 
Ck- 

In this article, our basic strategy is that we transform a 
complicated atomic configuration into PDs and try to identify 
meaningful shape information from specific distributions such 
as curves or islands in the PDs. Namely, we reconstruct char¬ 
acteristic atomic subsets from each distribution. To this aim, 
we compute the optimal cycle for each point (6^, dk) G Dg{A) 
on the distribution. Mathematically speaking, for a given ho¬ 
mology generator Ck of {bk,dk), the optimal cycle is obtained 
by solving a minimizing problem in the representatives of Ck 
under ^i-norm (see [18, 19]). Our method combining PDs 
with optimal cycles provides a tool to study inverse problems 
of PDs, and effectively works in the geometric analysis of glass 
structures, as we see shortly. 

For a mathematically rigorous introduction of these con¬ 
cepts, see Supporting Information or [10, II]. In this article, 
the PDs are computed by CGAL [20] and PH AT [21]. 


PDs for continuous random network structure 

Fig. 2 shows the PDs Di(Miiq), Di(Mamo), and Di(Ary) 
of a liquid Muq = (Qiiq, T^iiq), an amorphous Mamo = 
(Qamo ?-Ramo) 5 aud a Crystalline Mery — i^Qcry^Rcry) StatC of 


silica, respectively. Here, the horizontal and vertical axes are 
the birth and death scales, respectively, and the multiplicity of 
the PDs is plotted on a logarithmic scale. The configurations 
Qiiq, Qamo and Qcry arc acquired by molecular dynamics simu¬ 
lations using the BKS model [22]. The input radii R are set to 
be ro = 1.275 A and rsi = 0.375 A for each type of the atom 
(O or Si), which were determined from the first peak positions 
of the partial radial distribution functions of the amorphous 
configuration Qamo- The details of the molecular dynamics 
simulations and input radii are given in Supporting Informa¬ 
tion. 





a = ai 


a = a2 





Fig. 1. Atomic balls (top) and the PD Di(A) (bottom). New rings appear 
at ai{i = 2,3,4, 5), and the dashed rings express the disappearance. This is a 
schematic illustration showing the rings on Cp (red), Ct (blue), Cq (yellow), and 
Bq (green) in silica glass. The large and small balls correspond to oxygen and silicon 
atoms, respectively. 


We discovered that the PDs in Fig. 2 distinguish these three 
states. The liquid, amorphous, and crystalline states are char¬ 
acterized by planar (2-dim), curvilinear (I-dim), and island (0- 
dim) regions of the distributions, respectively. Here, the 0- and 
2-dimensionality of the PDs result from the periodic and ran¬ 
dom atomic configurations of the crystalline and liquid states, 
respectively. Furthermore, we emphasize that the presence of 
the curves in Di(Mamo) clearly distinguishes the amorphous 
state from the others. This implies that specific geometric 
features of the rings generating the curves in Di(Mamo) play 
a significant role for elucidating amorphous states. 

As shown in the figure, Di(Mamo) contains three character¬ 
istic curves Cp, Ct, and Co and one band region Ho, which 
are precisely characterized by using the invariance property 
with respect to the initial radius [15]. These particular dis¬ 
tributions, especially Co, start to become isolated near the 
glass transition temperature T = Tg (see Supporting Infor¬ 
mation) . Through further analysis of the persistent homology 
using optimal cycles, we found the following three geomet¬ 
ric characterizations, (i) The rings on Cp generate secondary 
rings on Ct,Co, and Bo (P is named after primary). That 
is, by increasing the parameter a, each ring on Cp becomes 
thicker and starts to create new rings by pinching itself, and 
these newly generated rings appear on Ct,Co, and Bq (a 
schematic illustration of the pinching process is described in 
Fig. I). This indicates a hierarchical structure from Cp to 
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Fig. 2. PDs of the liquid (left), amorphous (middle), and crystalline (right) states with the multiplicity on the logarithmic scale. In the amorphous state, the three 
characteristic curves and one band region are labeled Cp^CtiCq, and Bq, respectively. The insets in -Di(^amo) show rings in the hierarchical relationship, where the 
red and blue spheres represent oxygen and silicon atoms, respectively. 


Ct,Co, and Bq in the continuous random network. An ex¬ 
ample of the rings in this hierarchical relationship is depicted 
in the inset of Ili(^amo)- (h) The rings on Ct are constructed 
by tetrahedra consisting of four oxygen atoms at the vertices 
with one silicon atom at the center, (iii) The rings on Co 
and Bo are constructed only by the oxygen atoms (three and 
more, respectively). Recalling that the death scale indicates 
the size of rings, the rings on Ct are classified as SRO, while 
those on Cp, Co, and Bo are classified as MRO. 


Decomposition of FSDP 

The FSDP observed in the structure factor S{q) {q ~ 1.5-2 
A has been supposed to be a signature of MRO, but its 
real space origin is still controversial [4]. Here, we found that 
the distributions Cp, Co, and Bo of the MRO rings reproduce 
the q values of the FSDP fairly well. Moreover, we classified 
the MRO rings as a real space origin of the FSDP. 

We first note that the death scales of the rings on Cp, Co, 
and Bo are determined only by the oxygen atoms, and this is 
directly verified by the PD computation. In addition, recall 
that a is the parameter controlling the radius (a) = ar‘f 
of the z-th atomic ball Bi(a), and the death scale a — dk in¬ 
dicates the size of the individual ring Cfc. More precisely, the 
ring Ck disappears oX a — dk by being covered up in the atomic 
ball model [jf=i Bi{a). Hence £{dk) = 2y/dk + Tq measures 
the size of Ck on Cp, Co and Bo , where ro is the input radius 
of the oxygen atom. 

From the aforementioned argument, we define a distribution 


MaM) 
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i{dk)) 


where Ai = Cp U Co U Bo , A2 = Cp, A 3 = Co, and A 4 = Bo , 
and \Ai\ is the number of the elements in Ai. Here, S is the 
Dirac delta function which is used to count the contribution 
of each MRO ring in g-space. 

Fig. 3 shows the plots of MaA^) fhe structure factor 
S{q) around the FSDP. We found good agreement between 
the q values of the FSDP and the peak of MaA^). This im¬ 
plies that the MRO rings composed of Co, Cp, and Bo are 
the real space origin of the FSDP. We also note that the dis¬ 
tributions Mas{q)^ Ma 2 {q)^ MA 4 {q) are located on the large, 
medium, and small q values and, hence, the rings in Co, Cp, 
and Bo provide a decomposition of the FSDP into those q 
values, respectively. It should be emphasized that the MaA^) 
are derived from the configuration of the oxygen atoms only. 
This shows that the configuration of oxygen atoms plays a sig¬ 


nificant role in the FSDP. Furthermore, the invariant property 
[15] of i{dk) for the rings on Cp, Co and Bo induces that of 
MAi{q) under the choice of the input radius ro- This means 
that our analysis using PDs does not contain any artificial 
ambiguity of the input radii. 



Fig. 3. The distributions M^^iq) (z = 1,2, 3, 4) and the structure factor 
S{q) around the FSDP. 


Curves and shape constraints 

The presence of curves Cp, Ct, and Co in Di(Mamo) clearly 
distinguishes the amorphous state from the crystalline and liq¬ 
uid states. We emphasize here that this characteristic prop¬ 
erty shows the constraints on the shapes of the rings induced 
by the normal directions of these curves. 

For example, the shape of a ring on Co, which con¬ 
sists of three oxygen atoms (see the right panel in Fig. 4), 
is determined by specifying the first and second minimum 
edge lengths di and 0^2 (di < ^ 2 ) and the angle 0 between 
them, and hence is realized in a three dimensional parameter 
space. Then, the constraint for Co to be the curve requires 
that these three variables (di ,d 2 , 0 ) satisfy a certain relation 
/(di, ^ 2 , d) = 0, and hence provides a restriction on the shape 
of the 0-0-0 triangles. Fig. 4 shows a plot of (di,d 2 ,d) for 
the 0-0-0 triangles on Co in Di(Mamo), and we found a sur¬ 
face corresponding to this constraint /(di,d 2 ,d) = 0 in the 
parameter space. This demonstrates one of the medium-range 
geometric structures in the amorphous state. It is worth not¬ 
ing that this new characterization of MRO of 0-0-0 triangles 
in three dimensional parameter space cannot be derived by 
separately analyzing the conventional distributions of lengths 
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or angles, since each of them can only deal with the single 
parameter. 



Response under strain 

The presence of the curves also indicates variations in the 
shapes of the rings. That is, by following each curve along 
its tangential direction and studying its rings, we can observe 
the deformation of the rings. It is reasonable to suppose that 
these variations are due to thermal fluctuations, and hence 
the deformations of the ring configurations along the curves 
are probably softer than those in the normal directions. Con¬ 
sequently, the response under strain is expected to follow the 
same shape constraints. 


the original coordinates (blue), respectively. The contours of 
Co are depicted on the coordinates along the tangential di¬ 
rection. The bottom three panels show the projections of the 
histograms on their normal axes. 

The figures show that the contours of the strained state 
shift along the original curves Cp, Co, and Ct in T)i(^amo); 
i.e., downward in Cp, leftward in Co, and almost fixed in Ct, 
respectively. This is in contrast to a linear shrink, in which 
the contours simply move in the direction of decreasing both 
birth and death scales because of the uniform shrink of the sys¬ 
tem size. This strongly suggests that the configuration 
reflects the shape constraints and supports the expected me¬ 
chanical response of PDs. We also note from the figures of 
that the contour of Ct is almost fixed compared with 
those of Cp and Co- Since Cp and Co (Ct) are classified 
as MRO (SRO), this implies that MRO is mechanically softer 
than SRO. Here, we remark that we can observe a similar 
behavior of the PDs for a shear deformation. 

We have revealed that PDs encode information about elas¬ 
tic response, similar to how the radial distribution function 
encodes volume compressibility [23]. It should also be empha¬ 
sized that the curves or islands appear in the PDs of only the 
solids. This evidence suggests that these isolated distributions 
are related to the rigidity of the materials. This hypothesis 
follows from the fact that isolations represent geometric con¬ 
strains reflecting mechanical responses. Future numerical and 
theoretical studies to unravel this relationship would be of 
great value. 



0.15 
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"< 
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P 0.144 
0.142 

Fig. 5. Top: Contour plots of Cp, Co, and Ct for the original configuration 
Qamo of the amorphous state (black), the 1% strained state (red), and 

the artificial 1% linear shrink of the original coordinates (blue), respectively. 

Bottom: Projections of the histograms on the normal directions. 
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To verify this mechanical response of PDs, we performed 
simulations of isotropic compression for our amorphous state, 
and computed the PD of the strained state. The strain is set 
to be 1 % of the original volume, which is sufficiently small 
to satisfy a linear response relation with the stress. The three 
panels on the top in Fig. 5 show the contours of the histograms 
restricted on Cp, Co, and Ct for the original configuration 
Qamo of the amorphous state (black), for the strained state 
(red), and for the artificial 1 % linear shrink Qamo^’^ of 



Fig. 6. PDs for the Lennard-Jones system with the multiplicity on the logarithmic 
scale. The left and right panels correspond to Di and D 2 respectively. Top panels 
correspond to FCC crystal at T = 0.1 and bottom panels correspond to amorphous 
state at T = 10“^. 

PDs for random packing structures 

We next study the geometry of amorphous states close to ran¬ 
dom packing structures. In this case, we found that both 
Di and D 2 capture characteristic features of the amorphous 
structures. Fig. 6 shows the PDs Di(^cry) and Di(^amo) for 
i — 1,2 of the mono-disperse Lennard-Jones system in the 
crystalline and the amorphous states, respectively. The input 
radii r — ri = -- - = rwinR are set to be zero, since changing 
r only causes translations of the PDs for the single compo¬ 
nent system. The details of the simulation are explained in 
Supporting Information. 

Similar to the case of silica, the crystalline structure is char¬ 
acterized by the island distributions in the PDs (see top pan- 
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els in Fig. 6). They correspond to the regular triangles in 
^i(-^cry) and the regular octahedra, the regular tetrahedra 
and the quartoctahedra [24] in D 2 {Acry) of the FCC configu¬ 
ration. For the amorphous structure, the curves in Di(Aamo) 
and T)2(^amo) represent variations of triangles and tetrahedra, 
respectively. We also note that the isolation of the octahedral 
distribution is preserved well even in T)2(^amo)- Us peak is 
separated from the curve of the tetrahedra (see bottom right 
panel in Fig. 6), demonstrating the quantitative classification 
into two typical local structures of different-sized cavities. 

In random packings, it is known that the atomic configu¬ 
ration can be divided into dense packing regions built from 
tetrahedra and the complement that patches those regions to¬ 
gether [25]. In particular, the network structure of the dense 
packing regions characterizes the global connectivity beyond 
MRO, which has not yet been investigated in detail. Note 
that T)2(^amo) clearly separates the dense packing regions as 
the deformation curve of the tetrahedra and identifies the oc¬ 
tahedral island as the main component of the complement 
structure. From the Alexander duality in (e.g., [17]), the 
connectivity of the dense packing regions can be studied by 
the cavities of the complement. Namely, we first extract all 
octahedra from the octahedral island in T)2(Aamo) cind con¬ 
struct a new set of points Aoct — (Qoit? -Roct) by putting a 
point at the center of each octahedron (left of Fig. 7). Here, 
we set Roct to be zero, since the standard deviation of the 
octahedral sizes is very small. 

The right of Fig. 7 shows D2{Aoct)^ and we clearly observe 
that almost all cavities are located close to the diagonal. This 
means that the octahedral distribution rarely generates persis¬ 
tent cavities, and hence by the duality, we can conclude that 
the dense packing regions mostly construct a giant connected 
network. It should be remarked that the dual treatment is 
much easier and computationally more efficient than directly 
studying the intricate huge network structures. We also em¬ 
phasize that the analysis here of treating the octahedra as a 
new input is an iterative usage of the PD method, and can be 
an effective approach for studying multi-scale geometry. 


corresponding to octahedra appears in T)2(Aamo°^’^^°) and the 
characteristic curves are also observed in and 
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Fig. 8. PDs for CuZr alloys with the multiplicity on the logarithmic scale. The 
left and right panels correspond to Di and D 2 respectively. In the top panels, 
Di) ai'e described. The middle and bottom panels show the PDs of 
the atomic configurations of only Zr atoms in the alloys. The middle panels corre¬ 
spond to Cu 5 oZr 5 o alloy, and the bottom panels correspond to CuisZrgs alloy. 
The blue and green spheres represent copper and zirconium atoms, respectively. 



Fig. 7. Left: Red balls express the octahedra in the amorphous structure of the 
Lennard-Jones system, and the empty part corresponds to the dense packing regions. 
Right: PD D2{A^ct) of the left figure with the multiplicity on the logarithmic scale. 

As an example of multi-component systems, we also stud¬ 
ied metallic glasses composed of Cu and Zr [26], in particular, 
focusing on CusoZrso and CuisZrgs that display the different 
glass forming ability. The PDs for these alloys are shown in 
Fig. 8. The input radii for Cu and Zr are set to be 1.30 
A and 1.55 A, respectively, for the multi-component PDs (top 
panels in Fig. 8). These values are obtained by the same 
procedure as for the silica. Even in the multi-component sys¬ 
tem, the PDs basically show similar behaviors to those of the 
Lennard-Jones system. Specifically, the island distribution 



Fig. 9. The normalized pair-distance distributions (top panel) computed from 
and the radial distribution function for CusoZrso alloy (bottom 
panels) around the split second peak. The black line in the top panel was obtained 
by the pair-distances of the atoms in H C D 2 )- whereas the pink line 

corresponds to those in the complement of B. 

In the random packing structure, the split second peak of 
the radial distribution function has been supposed to be a 
signature of MRO [27]. The shaded region of the radial dis¬ 
tribution function in Fig 9 (bottom panel) shows the split 
second peak of CusoZrso- Then, we found that the pair- 
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distance distribution of the atoms of generators in B 
[0.92,1.05] X [1.76,2.01] C 7^2 ^-Iso shows a clear 

splitting of the distribution (black line in the top panel of 
Fig. 9) in the same length-scale. Here, B is chosen to be 
a region around the octahedral distribution. Meanwhile the 
pair-distance distribution of generators other than B shows a 
slight change there (pink line in Fig. 9). This means that 
the generators around the octahedral distribution play a sig¬ 
nificant role for the split second peak. Therefore, similar to 
the FSDP in the silica, this result demonstrates that the PDs 
classify the length-scale of MRO from other scales. 

We also studied the PDs of only the Zr component. The PD 
Di for Zr in CusoZrso (middle left panel in Fig. 8) represents 
the existence of a hierarchical MRO structure similar to that 
of the silica in Fig. 2, whereas the PD Di for Zr in CuisZrss 
(bottom left panel in Fig. 8) does not show any hierarchical 
curves. An example of the hierarchical rings in CusoZrso is 
depicted as the insets in the PD. 

We here remark that CusoZrso has higher glass forming 
ability than CuisZrss [28]. Relationships between the glass 
forming ability and the geometric structure are now actively 
studied (e.g., [29]). Then, this result suggests another pos¬ 
sibility that the presence of the hierarchical MRO structure 
extracted from PDs is also related to the glass forming ability 
of the alloy. In view of the results that the hierarchical MRO 
structure in PDs plays an important role for characterizing the 
glass states in Si02, this statement seems to be reasonable. In 
order to understand the geometric and topological formation 
of the glass, it will be a great challenge to clarify this new 
perspective. 

Conclusion 

We have presented that PD is a powerful tool for geomet¬ 
ric characterizations of various amorphous solids in the short, 
medium, and even further ranges. In this work, we have ad¬ 
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Supporting Information: Hierarchical 
structures of amorphous solids character¬ 
ized by persistent homology 


Computational homology 

We summarize the computational homological tools used in 
this article. For further mathematical and computational de¬ 
tails, please refer to [10]. 

Given a pair A = (Q,7^) of an atomic conhguration Q = 
(cci, X 2 ■ ■ • ,xn) and a set R = (ri, ^ 2 , • • •, 'Tat) of input radii, 
the persistence diagram captures hierarchical structures of 
the topological features such as rings and cavities in a one- 
parameter family of alpha shapes A(Q, R; a). Here, the alpha 
shape A(Q^R]a) is constructed as follows. Let = U^i 
be the Voronoi decomposition for Q, where Vi is the Voronoi 
cell of Xi. For each parameter a, we assign a ball Bi{a) = 
{cc G R^ I ||cc — cc^ll < ri(a)} with the radius ri{a) = 
for the zth atom. Then, the alpha shape A(Q, R; a) is defined 
by the dual of the cut balls Ci{a) = V D Bi{a), i— 1,..., 
Namely, this is a polyhedron with the vertex set Q such that 
an edge \xiXj \ (triangle, tetrahedron, respectively) is assigned 
to A{Q^R]a) if and only if the intersection Ci{a) f]Cj{a) of 
the corresponding double (triple, quadruple, respectively) cut 
balls is nonempty (see Fig. 10). In this article, the alpha 
shapes were computed using CGAL [20]. 

One of the important properties of the alpha shape is that 
the atomic ball model B{a) = alpha 

shape A{Q,R; a) can be continuously deformed to each other. 
Hence, there is no loss of topological information in the use 
of the alpha shape model, which is much easier to analyze 
computationally than the atomic ball model. Furthermore, 
we can trace the hierarchical structures in B{a) by changing 
the parameter a > amin, where amin = — min{ri, .. .r%}. In 
other words, the parameter a controls the resolution from fine 
(a ~ amin) to crude (a ^ <amin) geometric features. 

The left of Fig. 11 shows a schematic illustration of atomic 
ball models B{a) with the rings on Cp (red), Ct (blue), Co 
(yellow), and Bo (green). Here, the large (small) balls corre¬ 
spond to oxygen (silicon). At each ai {i = 2,3,4, 5), a new 
ring appears, and the dashed rings express the rings that have 
disappeared. The hierarchical ring structures are observed at 
0 ^ 3 , 0 ^ 4 , and as, at which new rings are generated by pinching 
the primary ring. 



The persistence diagram is a two-dimensional histogram 
recording the birth and death scales of the topological fea¬ 
tures in a one-parameter family {B(a)}a of atomic ball mod¬ 
els. For example, let c be the yellow ring in Fig. 11. We 
observe that c first appears a.t a = a 4 and first disappears 
at some value a*, where < a^. Then, the birth 

and death scales b and d, respectively, of the ring c are deter¬ 


mined by 6 = a 4 and d — a^. The histogram of the birth and 
death pairs (6, d) G R^ for all the rings appearing in the one- 
parameter family {B(a)}cx is called the persistence diagram 
Di{A). The right of Fig. 11 shows the persistence diagram 
of the left. From the definition of the persistence diagram, 
the birth scale indicates the maximum distance between adja¬ 
cent atoms in the ring, whereas the death scale indicates the 
size of the ring. It should also be mentioned that persistence 
diagrams Di(A) can be similarly defined for any nonnegative 
integer i G {0,1,2,...}, and they capture higher-dimensional 
topological features e.g., cavities for i — 2. The persistence 
diagrams in this article were computed by THAT [21]. 

MD simulation for the silica system 

The atomic configurations of silica in a liquid Quq, an amor¬ 
phous solid Qamo, and a crystalline solid Qcry state were gen¬ 
erated by the MD simulation. The system was composed of 
2,700 silicon atoms and 5,400 oxygen atoms interacting via the 
BKS potential [22], and the 24-6 Lennard-Jones potential is 
assigned to correct for the unphysical behavior during equili¬ 
bration in the liquid state [30]. The masses of the silicon and 
oxygen atoms were 27.98 and 15.99 grams/mole respectively. 
Despite its simplicity, the BKS model reasonably reproduces 
the static and dynamic properties of amorphous silica [31]. 

Starting from the beta-cristobalite structure, the crystalline 
configuration Qcry was obtained after the equilibration for 5 
ps at 1 bar and 10 K. The MD simulation was performed with 
Parrinello-Raman barostat and Nose-Hoover thermostat with 
the damping parameters Tbar = 0.5 ps and rth = 0.8 ps, re¬ 
spectively. The time step was set to be 1.0 fs. 

The liquid configuration Quq was obtained after the equi¬ 
librating the system from a random initial conhguration for 
50 ps at 1 bar and 7,000 K. Here, the MD simulation was 
performed with the Andersen barostat and the Nose-Hoover 
thermostat. 

The amorphous conhguration Qamo was obtained after the 
cooling MD simulation [31]. Starting from the liquid conhg¬ 
uration at the temperature 7,000 K, the heat-bath temper¬ 
ature was decreased linearly to the hnal temperature 10 K 
with —dTjdt = 10^^ K/s. Throughout the cooling process, 
the pressure was kept at 1 bar and we used the Nose-Hoover 
thermostat with rth = 1.0 ps. In the high temperatures re¬ 
gion 1,500 K < T < 6,000 K, we used Anderson barostat 
with Tbar = 1.0 ps. In the low temperatures region 10 K 
< T < 1,500 K, we used the Parrinello-Raman barostat with 
Tbar = 0.5 ps. Then, Qamo was obtained as the hnal conhg¬ 
uration of the cooling process. We tested hve cooling rates, 
—dT/dt = 10^^, 10^^, 10^^, 10^^, 10^^ K/s, and checked that 
these choices do not affect our results qualitatively. 

The set R = (ti,--- ,rw) of input radii can be adjusted 
for each zth atom (i = 1,...,^"). In our analysis, we chose 
the radii rsi and ro of silicon and oxygen using the hrst peak 
positions of the partial radial distribution functions of Qamo- 
The hrst peaks for SiO, GO, and SiSi appeared at rsio = 1.65 
A, Too = 2.55 A, and rsiSi = 3.15 A, respectively. Then, we 
selected the smallest two peaks and determined rsi and ro by 
solving rsi + To = Tsio and 2ro = too. This led to ro = 1.275 
A and rsi = 0.375 A. 

The glass transition temperature Tg and the melting tem¬ 
perature Tm are estimated through the enthalpy temperature 
curves. Starting from beta-cristobalite conhguration, the sys¬ 
tem is heated with a heating rate dT/dt = 0.1 K/ps from 
T = 1 K to T = 7, OOOK (black curve in Fig. 12). Then, Tm is 
evaluated to be around 3, OOOK [32]. Subsequently, the cooling 
simulation is performed with a cooling rate —dT/dt = 10.0 
K/ps (red curve). Then the glass transition temperature is 
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Fig. 11. Schematic illustration of atomic ball models showing the rings on Cp (red), (7 t (blue), Cq (yellow), and Bq (green). Here, the large (small) balls 

correspond to oxygen (silicon). As we increase the parameter oc, the primary ring generates the secondary rings in the hierarchy. 


evaluated as an intersection point between tangential lines of 
low (blue line) and high (green line) temperature regions of 
the curve and found to be Tg 2,700 K [31]. Fig. 13 shows 
PDs for 2, 500 K (left) and T = 3, 500 K (right), and the curve 
Co starts to appear below Tg. 



TK 

Fig. 12. Enthalpy as a function of temperature for cooling and heating of silica 
system. 
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Fig. 13. Di{A) below Tg (left) and in the super-cooled liquid region (right). 

MD simulation for the LJ system and the CuZr alloy 

The Lennard-Jones system is composed of 4,000 monatomic 
particles interacting via Lennard-Jones (LJ) potential u{r) — 
4£((cr/r)^^ — (cr/r)®) with £, a and mass equal to unity. The 
time step is set to be 0.005 in LJ units. 


Starting from the FCC configuration, is obtained after 
equilibration at temperature T = 0.1 and pressure p — 1 with 
Tth = 1.0 and Tbar = 1.0. To obtain an amorphous structure 
Qamoi fhe equilibration at T = 2 in the normal liquid phase 
has been performed for 2000 steps at the constant number 
density 1.015 [33] with Tth = 0.5. Then, the cooling MD sim¬ 
ulations have been performed for several cooling ratios from 
—dT/dt = 10“^ to 10“^. In the left panel in Fig. 14, the 
enthalpy normalized by the number of particles is described 
as a function of T. As is observed, there are no glass transi¬ 
tion in the LJ system. Therefore, we use an inherent structure 
for the liquid state as Qamo (black points in the left panel in 
Fig. 14). 



Fig. 14. Enthalpy as a function of temperature during cooling processes for 
Lennard-Jones system (left) and CuZr alloy (right), respectively. 

For the Cu-Zr alloy system, MD simulations have been per¬ 
formed using the embedded-atom method potential [34]. The 
masses of Cu and Zr atoms were set to be 63.54 and 91.22 
grams/mole, respectively. The number of particles was 16,000, 
in which there were 8,000 of both copper and zirconium atoms 
for CusoZrso, and 2,400 were copper atoms and 13,600 were 
zirconium atoms for CuisZrgs- 

In this case, we can observe the glass transition (the right 
panel in Fig. 14) during a cooling simulation with the quench 
rate 10^^ K/sec [26] down to T = 1 K keeping the pressure 
p = 0 atm after the equilibration at T = 17, 000 K. Then, 
we use the final configuration at T = 1 K as the amorphous 
configuration. 

All MD simulations in this article were performed using 
LAMMPS [35]. 






































